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Abstract 

A stochastic theory for a branching process in a neutron population 
with two energy levels is used to assess the applicability of the differ- 
ential self-interrogation Feynman-alpha method by numerically esti- 
mated reaction intensities from Monte Carlo simulations. More specif- 
ically, the variance to mean or Feynman-alpha formula is applied to 
investigate the appearing exponentials using the numerically obtained 
reaction intensities. 

1 Introduction 

Energy dependent aspects of neutron counting have been used to detect 
and determine the content of nuclear materials for several decades. The 
prime application is the Differential Die-Away Analysis (DDAA) method 
[H El El IU [5] . It is a deterministic method in which the time dependence of 
the detection rate of fast neutrons, after that an initial pulse of fast neutrons 
is injected to the sample, is used to determine the presence of fissile materials 
such as 235 JJ and 239 Pu. Its application is suitable when the fissile material 
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is embedded in a moderating surroundings such that the source neutrons 
induce thermal fission after having slowed down. 

The drawback of pulsed measurements such as the DDAA method is the 
necessity of the use of a neutron generator. The stochastic generalization of 
the DDAA method was recently suggested as an alternative method which 
eliminates the need for such an external source, namely the so-called Differ- 
ential Die-away Self-Interrogation (DDSI) technique [B]. The DDSI method 
utilizes the inherent spontaneous neutron emission of the sample. In the ab- 
sence of a trigger signal, the temporal decay of the correlations as a function 
of the time delay between two detections of fast neutrons is used in the DDSI 
method. This corresponds to a Rossi-alpha measurement with two energy 
groups. 

In Ref. [6J the functional form of the DDSI formula, i.e. the two-group 
Rossi-alpha formula, was derived in an empirical way. This formula was 
therefore later derived from first principles with the use of backward-type 
Kolmogorov equations (71 E]. A two-group theory of the Rossi- and Feynman 
alpha formulae is interesting also in areas other than nuclear safeguards, such 
as pulsed and stationary source driven experiments measuring the reactivity 
in fast cores of accelerator driven sub-critical systems. In such experiments 
it was found that two exponentials appear, indicating that the temporal 
behavior of the fast and thermal neutrons is separated, especially in fast 
reflected cores. This amplifies the need for the two group versions of the 
Feynman and Rossi-alpha formulae [5] . 

In this contribution, a stochastic theory for a branching process in a neu- 
tron population with two energy levels is investigated based on the previous 
results of Refs [7J El E] ■ In particular, a counterpart of the DDSI formula, the 
Differential Self-interrogation Variance to Mean or Feynman-alpha formula, 
which we shall call the DSVM formula, is derived by using the master equa- 
tion or Kolmogorov forward approach. The model includes a spontaneous 
fission source of fast neutrons, absorption in both groups, down-scattering 
(removal) from the fast to the thermal group, thermal fission, and detection 
of fast neutrons. 

Practical applicability of the method, which is based on the observation 
of two different exponentials, depends on the quantitative values of the var- 
ious within- and inter-group neutron reaction intensities, which appear as 
coefficients in the equations. This is in contrast to the energy-independent 
theory of multiplicity where the only appearing parameter is the first colli- 
sion probability. To assess the applicability and expected performance of the 
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DSVM method in practical situations, the above mentioned reaction inten- 
sities were determined from numerical Monte Carlo (MCNP4c) simulations. 
Significance of different values of the reaction intensities of thermal and fast 
neutrons in the performance of the method is discussed. 

2 A variance to mean formula for the de- 
tected fast neutrons 

In this section we will briefly discuss the fundamentals of the derivation of the 
variance to mean or Feynman-alpha for the two particle type system by using 
the Kolmogorov forward approach following Ref. [9]. In the model we have 
included a compound Poisson source of fast neutrons described by the source 
strength S\ that releases n particles with probability p q (n) at an emission 
event (i.e. spontaneous fission). We have assumed that the source is switched 
on at time t = t < 0, although the dependence on t will not be denoted. 
The detection rate of particles is included and is denoted by the intensity A^. 
We will start by giving a presentation of the analytical model consisting of 
a differential equation for the probability P(N±, N 2 , Z\, t) for having N\ fast, 
N 2 thermal neutrons at time t in the system and having detected Z\ fast 
particles in the interval (0,t). In deriving this differential equation we have 
summed all mutually exclusive events during an infinitesimally small time 
interval dt and we find for the probability P(N\, N 2 , Z\, t) the differential 
equation 



dt 



-(AiJVi + A 2 iV 2 + S 1 )P(N 1 , N 2 , Z u t) 



+ 



+ 



+ 



X la {N 1 + l)P{N 1 + l,N 2 ,Z 1 ,t) 
X 2a (N 2 + l)P(N 1 ,N 2 + l,Z 1 ,t) 
A^Aq + l^Aq + l,^-!,^,*) 



+ 



A 2/ (iV 2 + 1) J2 f(k)P(Ni -k,N 2 + 1, Z u t) 



k 



+ 



X d (N 1 + l)P(N 1 + l,N 2 ,Z 1 -l,t) 
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Here, Ai and A2 are the decay constants (total reaction intensities) for fast 
and thermal particles, whereas Ai a , A 2a are the absorption (actually, capture) 
intensities of fast and thermal particles, respectively. The removal of fast 
particles into the thermal group is described by Xr while fission resulting 
from the thermal particles happens with the intensity of A2/. The total 
intensities are given by 

Ai = Ai a + X R + A d , (2) 

and 

A2 = A 2a + A 2 /. (3) 

We derive the equations for the factorial moments by using the generating 
function of the form 

G(x,Y,z,t) = Y,J2z2 xNlYN2zZip ( N ^ N ^ z ^> ( 4 ) 

iVi N 2 Z x 

and describe the time evolution of the process by a partial differential equa- 
tion in terms of the generating function as, 

f)C f)C r)C 

— = (X la + X R Y + X d Z-X 1 X)— + (X 2a + X 2f q(X)-X 2 Y) — 

+ S 1 (r(X)-l)G, ' (5) 

where 

q(X) = ]T/ fc X fc , and (6) 

k 

r(X) = J>(n)jr\ (7) 

n 

Here, fk is the probability of having exactly k neutrons produced in an in- 
duced fission event. We have used the definition of the derivatives of the ex- 
pressions and ([7]) as Ui = dq/dX\ x =i and r x = dh/dX\ x =i, which stand 
for the expectations of the number of neutrons from induced and sponaneous 
fissions, respectively [TU]. We note that for t — > —00, the expectations of 
fast neutrons ((iVi)) and thermal neutrons ((N2)) will reach steady state due 
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to the stationary source term with intensity Si. The solutions to the sys- 
tem of differential equations are found by differentiation of equation (j5J) with 
respect to (X, Y, Z) and then letting (X — Y — Z — l). These read as 

A1A2 — VlARA2f (jJ\(jJ2 

(N 2 ) = N 2 = ^±, (9) 
(Zi) = eX 2f Nit, (10) 

where e = Xd/^2f and we have used the additional definitions Ui and w 2 and 

u eff as 

-wi = -^(Ai + A 2 ) + ^(Ai - A 2 ) 2 + 4A 1 A 2 ^ e// , (11) 

1„ .si 



-w 2 = -gCAi+Aa)- -^/(Ai-A 2 ) 2 + 4AiA 2 z/ e// , (12) 

^// = fi-r-jr"- (13) 

It can be noted that in reactor physics terminology, v e ff is identical with k e ff 
of a multiplying system and that the expectation of the detections increases 
linearly with time. This is due to the fact that it is determined by integration 
of the expectation of the neutron number with respect to time. In order to 
find the variance of the detector counts we need to determine the second 
moments by yet another differentiation with respect to (X, Y, Z) followed 
by letting (X = Y = Z = 1). We find the variance of the detector counts 
through the relation cr| = (Zi) + fx zz where the modified variance fizz is 
defined as fizz = (Z(Z — 1)) — (Z) 2 = a\ z — (Z) while in general we have 
Hxy = {XY) — (X)(Y). The differentiation procedure results in a system 
of six ordinary differential equations for the second order modified moments. 
The differentiation procedure gives a system of six dynamical equations of 
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the modified second moments as 

d 

^xx = -^Xxfixx + 2u 1 \ 2f HxY + ^X 2f N 2 + Str 2 , (14) 
d 

~q^ 1 xy = -(Ai + A 2 )/i X y + X R fi X x + fiA 2 /^yy, (15) 
.9 

— //yy = -2A 2 /iyy + 2X R fi X Y, (16) 
at 

T^ZX = -Xl^zX + VlX 2 fflZY + X d nxx, (17) 

T^/izY = -X 2 Hzy + X R /J,zx + X d HxY, (18) 

— nzz = 2e\ 2f fixz, (19) 

The equation system and its solution is rather analogous to the case of the 
Feynman-alpha equations with one neutron energy group but including de- 
layed neutrons, as given in Ref. [llj. Although an analytical solution for 
the general time-dependent system of equations ffT4l - f|T9|) would be hard 
to find, we note that in the stationary state the system breaks down into 
two systems such that the solution of the first three equations is indepen- 
dent from the second, such that the moments (f-ixx) = P>xx, (hxy) = fixY 
and (/iyy) = /iyy are constants. The equations describing detected particles 
need to be solved retaining the full time evolution by e.g. Laplace trans- 
forms. Moreover, it is found that the sought moment (fizz) is determined 
by quadrature of moment (nxz)- We find the constant 2nd moments as, 

(A| + UxU^v-iXyN-i + S x r 2 ) /nnS 

V>XX = ^7T ' 7 \ , (/UJ 

2(Ai + A2V1W2 
X 2 X R {v 2 X 2f N 2 + 5"ir 2 ) 

vxy - — oT\ — r~r^ ' \ zl ) 

2(Ai + A 2 )o;iC<; 2 
A|(^ 2 A 2/ iV 2 + S ir2 ) 

^YY = 7^77 ' 7 \ • \2l) 

2(Ai + A 2 )o;ia; 2 

Here we have used the notations v 2 = d 2 q/dX 2 \x=i and r 2 = d 2 h/dX 2 \x=i 
for the second factorial moments [10]. The objective now is to solve (fT7l) and 
([TBI by Laplace transform methods and we find the transformed identity as, 

viX d X 2 fjixY . (s + X 2 )X d jixx , no -> 

Hxz - 7j-T\ 1 77T\ { Z6 

sH(s) sH(s) 
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with 



H(s) = s 2 + (co 2 + wi)s + U!U 2 . (24) 

Note that we have assumed that the initial values of the moments \xxz and 
Hyz were equal to zero at t = (at the start of the measurement), hence 
the roots of H(s) determine the temporal behavior of the Feynman-alpha 
formula. Moreover, the solution has many similarities to that found in Ref. 
[TT] . We determine the variance to mean or Feynman-alpha formula by uti- 
lizing the relation for the variance azz = (Z) + fizz and after some algebra 
we find, 

w =1+yi(1 -^ )+y2(1 -^^ ) - (25) 

Here, the complete expressions for Y\ and Y 2 are quite lengthy. However, 
it turns out that the sum Yq = Y\ + Y 2 takes a rather simple form that 
also determines the value of the Feynman-alpha for large measurement times 
T — > oo as, 

Y = Y 1 + Y 2 = u 2 ^pl. (26) 

In the next section we will consider some quantitative examples of the vari- 
ance to mean formula in Equation f l25l) . 



3 Quantitative assessment of the variance to 
mean formula 

We will now discuss the quantitative application of the variance to mean for- 
mula (1251) for a MOX fuel assembly, considering setups when the MOX fuel 
was homogeneously mixed with a different amount of polyethylene modera- 
tor. The fuel sample was in form of a sphere which has the same amount of 
material as one single fuel pin, yielding a radius of 3.36 cm, of typical MOX 
fuel (U(94.14%)02 + Pu(5.86%)02) with approximately 4% fissile content. 
The amount of moderator mixed into the fuel was taken as that correspond- 
ing to an outer reflector shell of 4 different thicknesses: 5, 10, 15 and 20 
cm. In all cases the neutron source was assumed in the sample centre. In 
addition, a case was considered where the source was assumed to be equally 
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distributed in the volume (designated as VDS, Volume Distributed Source) 
with the amount of moderator corresponding to 5 cm thickness. The geom- 
etry of the calculations is displayed in the r.h.s. of Figure 1, whereas the 
l.h.s. of the Figure shows the equivalent heterogeneous model, displaying a 
surrounding moderator which was mixed homogeneously in the sample for 
the simulations. We have hence carried out five Monte-Carlo (MCNP4) sim- 
ulations [T2] in order to estimate the reaction intensities in the spent fuel 
assemblies with varying amount of moderation where the radius is taken 
as fissile material + moderator (3.36 + d cm). In the simulations we have 
assumed a two energy group approximation with a 0.46 eV energy cut-off. 
Moreover, the energies of the source neutrons were sampled from the Watt 
fission spectrum, originating from the spontaneous fission of 2m Pu. The 
energy spectrum is determined by the function, 

p(E) ~e- E/a smh(bEy/ 2 (27) 

where a = 0.799 MeV and b = 4.903 MeV- 1 . 



Table 1: Simulation Results 



Intensity 


5 cm 


10 cm 


15 cm 


20 cm 


VDS 


Ai [l/s] 


1.560 


1.230 


0.625 


0.185 


1.490 


A 2 [l/s] 


0.637 


1.220 


1.260 


1.180 


0.325 


\r [l/s] 


0.450 


0.621 


0.391 


0.094 


0.452 


A 2/ [l/s] 


0.243 


0.285 


0.171 


0.094 


0.118 



The simulation results are shown in Table [T]for the cases described above. 
The reaction intensities obtained from the simulations were normalized to one 
starting neutron and are obtained with a relative error less than 1%. The 
resulting reaction intensities are used in the variance-to-mean expression in 
Eq. ( 125]) and the results are displayed in Figures 2 A and B. In Figure 2 A, 
the variance to mean in lin-lin scale is shown for the parameters in Table 
1 whereas the remaining parameters are V\ = 2.87, v-i = 4.635, = 0.1, 
r% = 1.0 and r 2 = 0.0. We find that that the difference between the variance 
to mean for the volume distributed source (VDS in table ([1]), green line) and 
the point source with 5 cm moderator (blue line) is very small indicating 
that the point source model is sufficient for describing the system. We note 
that the values of v e ff are 0.316, 0.338, 0.244, 0.116 and 0.316 for the cases 
with 5 cm, 10 cm, 15 cm, 20 cm moderator and the VDS, respectively. In 
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POLYETHYLENE: CH 
Neutron source 
MOX fuel. 

d = (5-20) cm 



POLYETHYLENE: CH, + MOX fuel 




Heterogeneous case 



Homogeneous case 



Figure 1: Show the geometry used for the Monte Carlo (MCNP4c) simula- 
tions for the homogenized cases (left figure showing fuel mixed with different 
amounts of polyethylene moderator) compared to the heterogeneous cases 
(right). 



Figure A 



-5cm moderator 
-10cm moderator 
-15cm moderator 
-20cm moderator 
VDS 



-5cm moderator 
- 1 0cm moderator 
-15cm moderator 
-20cm moderator 
VDS 



100 
time 



Figure 2: (A and B) The Feynman- alpha formula is shown for the cases 
generated by Monte Carlo (MCNP4c) simulations with 5-20 cm moderator 
in lin-lin scale (Figure A) and lin-log scale (Figure B). 

Figure 2B the variance to mean from Figure 2A is displayed in lin-log scale 
using the same parameter values as in Figure 2A. The reaction intensities of 
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the moderated cases are approximately of the same order of magnitude and 
thus the two exponentials are not easily distinguished. 

4 Discussion and conclusions 

We have developed a forward Kolmogorov approach for the two group theory 
of the Differential Die- Away Self Interrogation Variance to Mean (DSVM), 
including a compound Poisson source and the detection process. The results 
agree with those calculated by the backward approach as reported in [7J [8] . 
We have used Monte Carlo simulations to find the reaction intensities needed 
to quantitatively assess the DSVM formula. We find that, unlike in the DDSI 
method (i.e. the two-group version of the Rossi-alpha method), the presence 
of two exponents in the solution is most often not clearly visible. This means 
that detection of the presence of fissile material may not be as obvious as with 
the Rossi-alpha method. By using the two group Rossi-alpha expression 4.34 
in Ref . [7J , we find the reverse situation, where in most cases with moderator 
the two exponentials appear. On the other hand, the determination of the 
exponents U\ and U2 by curve fitting could be more accurate in certain cases 
than with the DDSI method. Elucidating on the diagnostic value of the 
exponents in terms of determination of the sample parameters is not clear 
yet, and it requires further investigations, which will be reported in future 
work. 
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